% 2-D DIFFUSION-BASED RESPONSE MATRIX SOLVER
% TWO-DIMENSIONAL BENCHMARK PROBLEMS
% 1. IAEA 
% 2. BIBLIS
% 3. KOEBERG (4 group)
clear all
name = 'two_d_benchmarks';
prob = 2;

%prob    = 0;  % 1, 2, or 3 for IAEA BIBLIS or KOEBERG
locsolv = 0;  % 0 for KERNEL or 1 for FINITE DIFFERENCE
locfit  = 0;
order = 0; % legendre order
bc(1) = 1; % (L)
bc(2) = 0; % (R)
bc(3) = 1; % (B)
bc(4) = 0; % (T)
locsolv = 0;
fit=1;
in.reflector = 2;
if prob==1
%--------------------------------------------------------------------------
% PROBLEM 1: IAEA
%--------------------------------------------------------------------------
%--- MATERIALS
numtypes    = 4;   % number of element types
numg        = 2;   % number of groups
numm        = 4;   % number of materials
buck = 0.8e-4;
dat = [ 1.500   0.030   0.000  1.000  0.000  0.000  
        0.400   0.080   0.135  0.000  0.020  0.000
        1.500   0.030   0.000  1.000  0.000  0.000
        0.400   0.085   0.135  0.000  0.020  0.000
        1.500   0.030   0.000  1.000  0.000  0.000  
        0.400   0.130   0.135  0.000  0.020  0.000 
        2.000   0.040   0.000  0.000  0.000  0.000
        0.300   0.010   0.000  0.000  0.040  0.000 ];     
%dat(:,2)=dat(:,2)+buck*dat(:,1); % adjust for buckling   
abdat=0; % not used  
%--- PLACE THE ELEMENTS
%  elements = [ 3 2 2 2 2 2 2 3 3 2 2 2 2 1 1 4 4   % 1
%               2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 4 4   % 2
%               2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 4 4   % 3
%               2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 4 4   % 4
%               2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 4 4   % 5
%               2 2 2 2 2 2 2 2 2 2 2 1 1 4 4 4 4   % 6
%               2 2 2 2 2 2 2 2 2 2 2 1 1 4 4 4 4   % 7           
%               3 2 2 2 2 2 2 3 3 1 1 1 1 4 4 0 0   % 8
%               3 2 2 2 2 2 2 3 3 1 1 1 1 4 4 0 0   % 9
%               2 2 2 2 2 2 2 1 1 1 1 4 4 4 4 0 0   % 10
%               2 2 2 2 2 2 2 1 1 1 1 4 4 4 4 0 0   % 11
%               2 2 2 1 1 1 1 1 1 4 4 4 4 0 0 0 0   % 12
%               2 2 2 1 1 1 1 1 1 4 4 4 4 0 0 0 0   % 13       
%               1 1 1 1 1 4 4 4 4 4 4 0 0 0 0 0 0   % 14
%               1 1 1 1 1 4 4 4 4 4 4 0 0 0 0 0 0   % 15
%               4 4 4 4 4 4 4 0 0 0 0 0 0 0 0 0 0   % 16
%               4 4 4 4 4 4 4 0 0 0 0 0 0 0 0 0 0   % 17
%               ];    
%  elements = [ 3 2 2 2 2 2 2 3 3 2 2 2 2 1 1 4 4   % 1
%               2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 4 4   % 2
%               2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 4 4   % 3
%               2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 4 4   % 4
%               2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 4 4   % 5
%               2 2 2 2 2 2 2 2 2 2 2 1 1 4 4 4 4   % 6
%               2 2 2 2 2 2 2 2 2 2 2 1 1 4 4 4 4   % 7           
%               3 2 2 2 2 2 2 3 3 1 1 1 1 4 4 4 4   % 8
%               3 2 2 2 2 2 2 3 3 1 1 1 1 4 4 4 4   % 9
%               2 2 2 2 2 2 2 1 1 1 1 4 4 4 4 4 4   % 10
%               2 2 2 2 2 2 2 1 1 1 1 4 4 4 4 4 4   % 11
%               2 2 2 1 1 1 1 1 1 4 4 4 4 4 4 4 4   % 12
%               2 2 2 1 1 1 1 1 1 4 4 4 4 4 4 4 4   % 13       
%               1 1 1 1 1 4 4 4 4 4 4 4 4 4 4 4 4   % 14
%               1 1 1 1 1 4 4 4 4 4 4 4 4 4 4 4 4   % 15
%               4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4   % 16
%               4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4   % 17
%               ];              
% elements = [3 2 2 2 3 2 2 1 4
%             2 2 2 2 2 2 2 1 4
%             2 2 2 2 2 2 1 1 4
%             2 2 2 2 2 2 1 4 4
%             3 2 2 2 3 1 1 4 4
%             2 2 2 2 1 1 4 4 4 
%             2 2 1 1 1 4 4 4 4
%             1 1 1 4 4 4 4 4 4 
%             4 4 4 4 4 4 4 4 4
%           ];
elements = [3 2 2 2 3 2 2 1 3
            2 2 2 2 2 2 2 1 3
            2 2 2 2 2 2 1 1 3
            2 2 2 2 2 2 1 3 3
            3 2 2 2 3 1 1 3 3
            2 2 2 2 1 1 3 3 3 
            2 2 1 1 1 3 3 3 3
            1 1 1 3 3 3 3 3 3 
            3 3 3 3 3 3 3 3 3
          ];      
      %elements=elements*0.0 + 1;
xcm  = [   0     20  ];
xfm  = [     1       ]*100;  
ycm  = [   0     20  ];
yfm  = [     1       ]*100;
% compute number in each region
numycm = length(yfm); % number of y regions
numxcm = length(xfm); % number of x regions
mt = zeros(numxcm,numycm,numtypes);
% ELEMENT 1 - material specification
mt(:,:,1)   = [  1 ]; 
% ELEMENT 2 - material specification
mt(:,:,2)   = [  2  ]; 
% ELEMENT 3 - material specification
mt(:,:,3)   = [  3 ]; 
% ELEMENT 4 - material specification
mt(:,:,4)   = [  4  ]; 

elseif prob ==2
    
%--------------------------------------------------------------------------
% PROBLEM 2: Biblis
%--------------------------------------------------------------------------
%--- MATERIALS
numtypes    = 8;   % number of element types
numg        = 2;   % number of groups
numm        = 8;   % number of materials
dat = [ 1.4360000  0.0272582  0.0058708  1.0000000  0.0000000  0.0000000
        0.3635000  0.0750580  0.0960670  0.0000000  0.0177540  0.0000000
        1.4366000  0.0272995  0.0061908  1.0000000  0.0000000  0.0000000
        0.3636000  0.0784360  0.1035800  0.0000000  0.0176210  0.0000000
        1.3200000  0.0257622  0.0000000  0.0000000  0.0000000  0.0000000
        0.2772000  0.0715960  0.0000000  0.0000000  0.0231060  0.0000000
        1.4389000  0.0274640  0.0074527  1.0000000  0.0000000  0.0000000
        0.3638000  0.0914080  0.1323600  0.0000000  0.0171010  0.0000000
        1.4381000  0.0272930  0.0061908  1.0000000  0.0000000  0.0000000
        0.3665000  0.0848280  0.1035800  0.0000000  0.0172900  0.0000000
        1.4385000  0.0273240  0.0064285  1.0000000  0.0000000  0.0000000
        0.3665000  0.0873140  0.1091100  0.0000000  0.0171920  0.0000000
        1.4389000  0.0272900  0.0061908  1.0000000  0.0000000  0.0000000
        0.3679000  0.0880240  0.1035800  0.0000000  0.0171250  0.0000000
        1.4393000  0.0273210  0.0064285  1.0000000  0.0000000  0.0000000
        0.3680000  0.0905100  0.1091100  0.0000000  0.0170270  0.0000000];
abdat=0; % not used  
%--- PLACE THE ELEMENTS   
elements = [ 1 8 8 2 2 6 6 1 1 7 7 1 1 4 4 3 3   % 1
             8 1 1 8 8 2 2 8 8 1 1 1 1 4 4 3 3   % 2
             8 1 1 8 8 2 2 8 8 1 1 1 1 4 4 3 3   % 3
             2 8 8 1 1 8 8 2 2 7 7 1 1 4 4 3 3   % 4
             2 8 8 1 1 8 8 2 2 7 7 1 1 4 4 3 3   % 5
             6 2 2 8 8 2 2 8 8 1 1 8 8 4 4 3 3   % 6
             6 2 2 8 8 2 2 8 8 1 1 8 8 4 4 3 3   % 7           
             1 8 8 2 2 8 8 2 2 5 5 4 4 3 3 3 3   % 8
             1 8 8 2 2 8 8 2 2 5 5 4 4 3 3 3 3   % 9
             7 1 1 7 7 1 1 5 5 4 4 4 4 3 3 3 3   % 10
             7 1 1 7 7 1 1 5 5 4 4 4 4 3 3 3 3   % 11
             1 1 1 1 1 8 8 4 4 4 4 3 3 3 3 3 3   % 12
             1 1 1 1 1 8 8 4 4 4 4 3 3 3 3 3 3   % 13       
             4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 3 3   % 14
             4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 3 3   % 15
             3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3   % 16
             3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3   % 17
             ];    
xcm  = [   0     20.1226/2  ]*2;
xfm  = [     1       ]*100;
ycm  = [   0     20.1226/2  ]*2;
yfm  = [     1       ]*100;
% compute number in each region
numycm = length(yfm); % number of y regions
numxcm = length(xfm); % number of x regions
mt = zeros(numxcm,numycm,numtypes);
% ELEMENT 1 - material specification
mt(:,:,1)   = [  1 ]; 
% ELEMENT 2 - material specification
mt(:,:,2)   = [  2  ]; 
% ELEMENT 3 - material specification
mt(:,:,3)   = [  3 ]; 
% ELEMENT 4 - material specification
mt(:,:,4)   = [  4  ];    
% ELEMENT 5 - material specification
mt(:,:,5)   = [  5 ]; 
% ELEMENT 6 - material specification
mt(:,:,6)   = [  6  ]; 
% ELEMENT 7 - material specification
mt(:,:,7)   = [  7 ]; 
% ELEMENT 8 - material specification
mt(:,:,8)   = [  8  ];    

elseif prob == 3
    
%--------------------------------------------------------------------------
% PROBLEM 3: Koeberg
%--------------------------------------------------------------------------
%--- MATERIALS
numtypes    = 7;   % number of element types
numg        = 4;   % number of groups
numm        = 7;   % number of materials
dat = [ ...
    2.491869  0.067927  0.008228  0.745248  0.000000  0.000000  0.000000  0.000000
    1.045224  0.066508  0.000536  0.254328  0.063789  0.000000  0.000000  0.000000
    0.677407  0.070757  0.007058  0.000424  0.000486  0.064381  0.000000  0.001245
    0.375191  0.069235  0.083930  0.000000  0.000000  0.000003  0.050849  0.000000
    2.492653  0.067275  0.008295  0.745248  0.000000  0.000000  0.000000  0.000000
    1.049844  0.065296  0.000713  0.254328  0.063112  0.000000  0.000000  0.000000
    0.676610  0.070432  0.009230  0.000424  0.000478  0.063078  0.000000  0.001543
    0.379481  0.086595  0.108244  0.000000  0.000000  0.000003  0.048420  0.000000
    2.491978  0.066922  0.008285  0.745248  0.000000  0.000000  0.000000  0.000000
    1.051910  0.064628  0.000713  0.254328  0.062765  0.000000  0.000000  0.000000
    0.677084  0.069952  0.009214  0.000424  0.000473  0.062404  0.000000  0.001598
    0.381453  0.089675  0.108087  0.000000  0.000000  0.000003  0.047549  0.000000
    2.492535  0.067963  0.008459  0.745248  0.000000  0.000000  0.000000  0.000000
    1.045298  0.066632  0.000923  0.254328  0.063737  0.000000  0.000000  0.000000
    0.674684  0.072139  0.011714  0.000424  0.000486  0.064330  0.000000  0.001630
    0.374240  0.092630  0.133600  0.000000  0.000000  0.000003  0.049518  0.000000
    2.492329  0.066940  0.008409  0.745248  0.000000  0.000000  0.000000  0.000000
    1.051953  0.064694  0.000921  0.254328  0.062737  0.000000  0.000000  0.000000
    0.675683  0.070681  0.011675  0.000424  0.000473  0.062376  0.000000  0.001797
    0.380606  0.102043  0.134282  0.000000  0.000000  0.000003  0.046859  0.000000
    2.491521  0.066584  0.008400  0.745248  0.000000  0.000000  0.000000  0.000000
    1.054029  0.064020  0.000921  0.254328  0.062386  0.000000  0.000000  0.000000
    0.676197  0.070201  0.011651  0.000424  0.000468  0.061696  0.000000  0.001852
    0.382438  0.105135  0.133974  0.000000  0.000000  0.000003  0.046005  0.000000
    2.119737  0.042840  0.000000  0.745248  0.000000  0.000000  0.000000  0.000000
    0.980098  0.044852  0.000000  0.254328  0.042052  0.000000  0.000000  0.000000
    0.531336  0.056528  0.000000  0.000424  0.000322  0.044589  0.000000  0.000978
    1.058029  0.117896  0.000000  0.000000  0.000000  0.000000  0.052246  0.000000
];
abdat=0; % not used  
%--- PLACE THE ELEMENTS   
elements = [ 1 3 3 1 1 3 3 1 1 2 2 1 1 4 4 7 7   % 1
             3 1 1 3 3 1 1 2 2 1 1 6 6 4 4 7 7   % 2
             3 1 1 3 3 1 1 2 2 1 1 6 6 4 4 7 7   % 3
             1 3 3 1 1 2 2 1 1 3 3 4 4 7 7 7 7   % 4
             1 3 3 1 1 2 2 1 1 3 3 4 4 7 7 7 7   % 5
             3 1 1 2 2 1 1 3 3 5 5 4 4 7 7 0 0   % 6
             3 1 1 2 2 1 1 3 3 5 5 4 4 7 7 0 0   % 7           
             1 2 2 1 1 3 3 1 1 4 4 7 7 7 7 0 0   % 8
             1 2 2 1 1 3 3 1 1 4 4 7 7 7 7 0 0   % 9
             2 1 1 3 3 5 5 4 4 7 7 7 7 0 0 0 0   % 10
             2 1 1 3 3 5 5 4 4 7 7 7 7 0 0 0 0   % 11
             1 6 6 4 4 4 4 7 7 7 7 0 0 0 0 0 0   % 12
             1 6 6 4 4 4 4 7 7 7 7 0 0 0 0 0 0   % 13       
             4 4 4 7 7 7 7 7 7 0 0 0 0 0 0 0 0   % 14
             4 4 4 7 7 7 7 7 7 0 0 0 0 0 0 0 0   % 15
             7 7 7 7 7 7 7 0 0 0 0 0 0 0 0 0 0   % 16
             7 7 7 7 7 7 7 0 0 0 0 0 0 0 0 0 0   % 17
             ];    
xcm  = [   0     21.608/2  ];
xfm  = [     1       ]*10;
ycm  = [   0     21.608/2  ];
yfm  = [     1       ]*10;
% compute number in each region
numycm = length(yfm); % number of y regions
numxcm = length(xfm); % number of x regions
mt = zeros(numxcm,numycm,numtypes);
% ELEMENT 1 - material specification
mt(:,:,1)   = [  1 ]; 
% ELEMENT 2 - material specification
mt(:,:,2)   = [  2  ]; 
% ELEMENT 3 - material specification
mt(:,:,3)   = [  3 ]; 
% ELEMENT 4 - material specification
mt(:,:,4)   = [  4  ];    
% ELEMENT 5 - material specification
mt(:,:,5)   = [  5 ]; 
% ELEMENT 6 - material specification
mt(:,:,6)   = [  6  ]; 
% ELEMENT 7 - material specification
mt(:,:,7)   = [  7 ]; 
% ELEMENT 8 - material specification
mt(:,:,8)   = [  8  ];    
    
end

tol     = [1.e-9 1.e-9]; % tolerance for norms             
maxit   = 20;

swpi = 1; % currently not used
swnt = 1; % currently not used
swgm = 0; % ''

input = makeinput( name, numtypes, numg, numm, order, dat, abdat, ...
                   elements, bc, xcm, xfm, ycm, yfm, mt, tol, maxit, ... 
                   swpi, swnt, swgm, 1 ) ; 
%input.x0 = x1;
input.locsolv = locsolv;
input.locfit = locfit;
%-------------------------------------------------------------------------
% RUN SOLVER and OUTPUT EITHER TO FILE OR SCREEN      
input.reflector = 2;

input.locsolv = 0;
input.fit = fit;
[x1,normk1,it,nrm1,fr] = solverPI_mc(input);
disp(['  iteration number: ', num2str(it), ' norm(x) = ', num2str(nrm1(it))])
keff=x1(end-1);
fprintf(1,'  final keff = %12.8f\n',x1(end-1));
if (locsolv==0)
    km=keff;
    frm=fr;
else
    kr=keff;
    frr=fr;
end


% rel=(efr0-efr0f)./efr0; rel=rel(rel<1000); 
% aa(ii) = (k0-k0f)*1e5; bb(ii) = max(abs(rel)); cc(ii) = mean(sqrt(rel.^2));
% fprintf(' %12.10f  %12.10f  %12.10f \n', aa(ii),bb(ii),cc(ii) )

%toc
% 
% tic
% input.maxit=1;
% [xseed,normk2o,it2o,nrm2o]=solverPI_mc(input);
% [x2,normk2,it2,nrm2] = solverNT_mc(input,xseed);
% disp(['  iteration number: ', num2str(it2), ' norm(x) = ', num2str(nrm2(it2))])
% fprintf(1,'  final keff = %12.8f\n',x2(end-1));
% toc
% it2  = it2o+it2;
% nrm2 = [nrm2o nrm2];
% semilogy(1:it,nrm1,'-o',1:it2,nrm2,'--x','LineWidth',2)
% legend('PI','Newton')
% xlabel('(outer) iteration')
% ylabel(' ||f(x)||_2 ')
% title('comparison of convergence rates')

